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Abstract 

Based on our previous work for solving the nonlinear Schrodinger equation with multichannel dynamics that 
is given by a localized standing wave and radiation, in this work we deal with the multichannel solution 
which consists of a moving soliton and radiation. We apply the modulation theory to give a system of 
ODEs coupled to the radiation term for describing the solution, which is valid for all times. The modulation 
equations are solved accurately by the proposed numerical method. The soliton and radiation are captured 
separately in the computation, and they are solved on the translated domain that is moving with them. 
Thus for a fixed finite physical domain in the lab frame, the multichannel solution can pass through the 
boundary naturally, which can not be done by imposing any existing boundary conditions. We comment on 
the differences of this method from the collective coordinates. 

Keywords: moving soliton, radiation, multichannel dynamics, nonlinear Schrodinger equation, modulation 
equations, numerical method, boundary condition 


1. Introduction 

We consider the following nonlinear Schrodinger (NLS) equation in d dimensions (d = 1,2,3) 

id t u(x, t) + Au(x, t) + (3(\u(x, t)\ 2 )u(x, t) = 0, x e M d , t > 0, (1) 

u(x, 0) = izo(x), x G M d , (2) 

where uo(x) is the given initial data and /?(•) : M -+ M is a smooth nonlinear function. It is well-known that 
the mass of the system M and the Hamiltonian (energy) of the system H are conserved, i.e. 


M(t) := 

f | u(x, t) 2 dx = M(0), 

( 3 ) 


J R d 


m ■■= 

[ 1" Vrt(x,£) 2 -F (V(x,£)| 2 )] dx = H(0), 

J Rd L \ yj 

( 4 ) 


with F'(p) = /3(p). We study the multichannel dynamics in §. In fact, many physical systems such as the 
particles or matter wave dynamics in the quantum mechanics and the nonlinear optics [32, 2, 25, 24, 26, 31 , 
involve the dynamics of the multichannel solution which mean the solution of the system asymptotically 
is given by a linear combination of a localized (in space), periodic (in time) wave (solitary or standing 
wave) and a dispersive part [35] . The dispersive part is usually referred as radiation in the literature. The 
multichannel solutions also widely exist in many other conservative nonlinear dispersive and wave equations 
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besides the NLS equation [35], and they have been found of importance and useful in both theoretical 
analysis and applications. 

In this paper, we focus on the study of the multichannel solution of a moving soliton and radiation in 
the NLS equation 0 . ie - the solution u(x, t ) of 0 is given by 

^(X, t) = $ soliton (X, t) + <f> radiation (x, t) , (5) 

where radiation^ t) denotes the dispersive wave and the soliton solitonfc t) has form 

$soHton(x, t) = e* 0(x>t) 0 w (x — vt — D), X G M d , t > 0, 

with vGl d the velocity, DGM d the shift in space, 0(x, t) £ M the phase and <p u (x) > 0 the eigenfunction 
of cj > 0 , i.e. 

W0 w (x) = A0 w (x) + ^(^(x))0 w (x), X e 
When the system is free from radiation, the phase function is given by 

„/ v v • x Iv1 2 t . . 

d ( x ,t) = — - — +w£ + 7, (6) 

where 7 G M is the shift in phase. In the presence of radiation, it then becomes a function in general and 
the parameters v,D,u; become functions of time. 

In the studies of the multichannel dynamics, a very popular method in the physics community is the 
collective coordinates method. The collective coordinates method usually begins with a guess of the soliton 
of the system and drops the radiation, which will result in some ODEs to describe the physical system 
approximately. Thus it often applies when the soliton part is the main interest and the radiation is small. 
For a detailed review of this method, we refer the readers to |31] and the references therein. However as 
pointed out in [39], this method can only provide a good approximation of the soliton for short time. For 
the long time dynamics or large initial radiation, the collective coordinates method will fail. Another major 
approach to study the multichannel dynamics is to numerically solving the governing equations by truncating 
the whole space problem onto a finite domain and then imposing some suitable boundary conditions. The 
imposed boundary condition for computation is the key to maintain the true physics. Classical Dirichlet 
or Neumann boundaries introduce spurious reflections, while periodic boundaries allow outgoing waves 
to wrap around the computational domain. For the wave equations, the Dirichlet-to-Neumann method 
that attempts to use the exact solution as a boundary condition, has been established and works very 
well g Emma mi nu. However it is fraught with problems for dispersive equations na im ed and 
only limited progress has been made so far. People then used a dissipative term which is localized on a 
buffer region to dissipates outgoing waves, but this method also dissipates incoming waves located near the 
boundary which is spurious HE]. By decomposing the solution into a family of coherent states, a phase 
space filter method has been proposed to design an open boundary for the NLS equation in [34, 33 , but 
this method fails to filter waves with wavelength longer than the buffer region, which is also the problem 
shared with most absorbing boundary conditions. In |33j, by combining the phase space filter method to a 
spectral technique that resolves waves of both long wavelength and short wavelength, a multiscale method 
has been considered to filter outgoing waves regardless of the frequency. 

To study the multichannel solution © in 0 , we are going to apply the modulation theory, which has 
been established in m for studying the stabilization of solution in the NLS equation from the asymptotical 
point of view. Here the modulation equations for governing the multichannel solution are derived exactly. 
They are a system of ODEs coupled to a dispersive equation for the radiation and they can exactly describe 
the dynamics of the solitary wave and the dispersive part separately at the same time, which are valid for 
all times. Thanks to our recent work in [39] where we studied the multichannel dynamics of a standing wave 
and radiation, an efficient and accurate numerical method is proposed to solve the modulation equations. 
The modulation equations are indeed solved numerically on the Lagrangian domain that is moving with 
the soliton. Thus, for any fixed finite domain in the lab frame, the multichannel solution can approach and 
pass through the boundaries naturally, which can not be done by any boundary techniques as we discussed. 
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Numerical explorations are provided in the end towards better understanding of the multichannel dynamics. 
Although the work here is done for the NLS equation, it is believed that the approach could work for other 
dispersive equations. 

The rest of the paper is organized as follows. In Section [2j we shall derive the modulation equations 
with a brief review. In Section [3J we shall propose an algorithm for the modulation equations. Numerical 
studies and results are given in Section [4] followed by some discussions. 


2. Modulation equations 

In this section, we shall formally derive the full modulation equations of the NLS equation for describing 
the multichannel solution 0 by following the sketch used in [ 12 ] and review some related mathematical 
theories for the readers’ convenience. Comments on the differences of this method from the well-known 
method of collective coordinates are given in the end. 

2.1. Formal derivation 

Based on the fact of the free soliton ([ 6 | in the NLS equation ([l]) and the physical observation that the 
qualitative behavior of the linear Schrodinger equation should not change that much in response to a small 
nonlinear and Hamiltonian perturbation in the dynamics, i.e. we should still see a localized part which 
decouples after a long time from the dispersive part, thus we take the ansatz of the solution of § m as 

w(x, t) = e^ (x>t) [0 w(t) (y(x, t)) + R (y(x, t),t)] , x e R d , t > 0, (Al) 

0(x,i) = ^v-x -\J l v (s)| 2 ds + ^ w(s)ds + 7 (i), (A2) 

y(x,i) = x - [ v(s)ds - D(t), (A3) 

Jo 

with initial conditions 

W 0 (x) = e*(~V +7 °) [<^ o (x - Do) + i? 0 (x - Do)], 

v(0) = v 0 , w(0) = w 0 , 7 ( 0 ) = 70 , D( 0) = D 0 , i?(x, 0) = i? 0 (x). (7) 

Here, 0 a; (x) is the nonlinear bound state of the time-independent NLS equation of 0 with eigenvalue cj, 

i.e. 

W0 w (x) = A0 w (x) + /3(^(x))0 w (x), X e (8) 

<Mx) G H\R d ), </> u > 0. 

y G is interpreted as the translated space variable or the Lagrangian domain which is moving with 
the solution with velocity v. 0 w (y) is the moving soliton in the multichannel solution and R(y,t) is the 
dispersive wave. Under some conditions on the nonlinearity /?(•), for uj > 0, the eigenvalue problem § has 
a unique radially symmetric solution exponential decaying at far field [HISS]. 

Plugging the ansatz into the NLS equation ([!]), we obtain 

id t R + ^yy — 2 .D • R — + y H—^ R — it) • d y 0^, + iivd^fp^ (9) 

+ 9yy(pu) — + 7 + ^ 0a; + 0 (10a; + A!| 2 ) (0a; + R) = 0, X G M^, t > 0, 

where R = i2(y(x, t), t), 0^ = 0 w ( t )(y(x, t)), d u <)) u = 9 a; 0 a; ( t )(y(x, t)), d y and d yy denote the gradient and 

Laplacian operators with respect to the variable y respectively, and the function 9 u; 0 u; (x) is given by 

ud w 0 w (x) = [A + 20'(0^(x))0 2 (x) + 0(0^ (x))] d w 0 w (x) - 0 w (x), x e M d . (10) 
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Since ^(x) is even, so it is clear that ^^(x) is also an even function. By introducing a new parameter 
£(t) with arbitrary initial value £(0) = £o G M and 


£(*) = i(t) + 2 




v(s)ds + -v(t) • D(t), 


t > 0, 


( 11 ) 


which indicates that 




y-y 

2 ’ 


the explicit dependence on the variable x in ^ can be removed, and the equation can be interpreted as 
imposed on the y-domain, i.e. the translated Lagrangian domain in terms of y. By further decompose the 
interaction term into the linear part of R and the nonlinear part denoted as N = N( y, t), i.e. 


P (l<^ + R\ 2 ) {K + R)=P {Pi) P» + P {Pi) R + P' {Pi) Pl-{R + R) + N, (12) 


where here and after z denotes the complex conjugate of a complex number z. The nonlinear part N = 
0(\R\ 2 ), as \R\ —y 0. For the cubic nonlinearity case, i.e. /3(p) = Ap, A G M, we have 

N = \<j> u {R 2 + |i?| 2 ) + \(p u + R)\R\ 2 . 

Then by also noting (J8|, the equation can be rewritten as 

id t R(y, t ) + (d yy - it) ■ <9 y ) R(y, t) - R{ y, *) (13) 

+ [P{Pl{ y)) + P'{Pl{y))Pl{y)\ R{y,t) + P'{Pl{y))Pl{y)R{y,t) + N( y , t) 

- it) ■ dypupy) + iud u (j> u (y) ~ (i + 0w(y) =0, y eR d , t > 0. 


As a convention [12, 30 ], in the following, we interpret R = Ri +iR-> G C which is originally a complex-valued 
scalar function, as a two-dimensional column vector function and so do other terms in (|13|) , i.e. 


R{ y,t) 


( Ri{y,t)\ 

\R2(y,t)J ’ 


PUy) 



N(y,t) 


f Ni{y,t)\ 

VlV 2 (y ,t)J’ 


where Ri and R 2 denote the real part and imaginary part of R respectively and Ni and N 2 denote the real 
part and imaginary part of N respectively. Introduce two two-by-two matrixes 


J = 



m = 


0 L- 


with time dependent operators 


L+ ■— ~dyy + W - P(pl ) - ZPPpl)p\ 


z, 


£_ := 


-8yy + U ~ P(PI). 


Here the matrix J satisfying J 2 = —Id is introduced as the matrix representation of the imaginary unit —i. 
Then (13) can be written into a real-valued vector form as 


9tR{y, t)=JH(t)R{ y , t) + {t> ■ dy j R(y , t) + 

+ Jpuii y) + (t> ■ dy j 


JR(y, t) - JN(y, t) 
Puj(y) - ud u pu{y), y e R d , t>o. 


(14) 


To make ( fl4| ) solvable, we need to impose some extra conditions on the soliton and the radiation. We assume 
R G {4>u, Jdu4>u, y4>uj, Jdy4 > u>}~ L , which is known as the orthogonality conditions in the literature mm 
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l22l 143] . and here the inner product of two vectors are interpreted as usual, 
conditions read 


(Ri{t, ■), = o, 

(Ri{t,-), y = 0 , 


(R-2(t,-), d u (j> u ( t )) = 0, 
(R 2 (t,-), <9 y <?W)) = 0 , 


In details, the orthogonality 

t > 0, (15) 


where < f,g >:= f Rd /(x) • g(x)dx for two real-valued scalar or vector functions /,g G L 2 (R d ). 

Now we firstly take the inner product of (14) with <p UJ ( y) on both sides. By noting that < 0 W , (D 
d y )0 W >= £>• < 0^, d y 0^ >= 0, we get 


^ (0w - #1, du><l>i 




yi?2) + D • (0^, 9yi?l) + £ (0u;> ^ 2 ) — (0ui5 ^ 2 ) 5 t > 0. 


(16) 


Secondly, noting the fact that 


(<9c>u,, L+iii) = - ( 0 W , i?i) = 0, 


which is indicated by integration by parts and (10), then by taking the inner product of (14) with Jd u (j) u { y) 
on both sides, we get 

i {{d u (f> u , Ri) + {4>u, d u <j>J)) =W (d^<j) u , R 2 ) - ^ • {{d u <j) u , yRi) + y^ w )) 

+1) • (d^, d y R 2 ) + (d^, JVi), t > 0. (17) 

Here <9^0^ is given by 

[w - A - 2f3'(4>l)(t>l - P(4>1)\ dl4> u = [-2 + 4 P"(4>l)4>ld u 4> u + 6/?'(<^)<£ w d w &,] d u <f> u . (18) 

Thirdly, by taking the inner product of ( |l4| with Jd y (f) u ,(y) on both sides and noting that < d y ^> u ,y >= 
~^\\(f>u\\ 2 L 2 , we get 

v I l (dy^, yRi) ~ lll^llls) = — to (dufa, d y R 2 ) + D {dycpu, d y R 2 ) - £ (<9 y <^, i?i) 

— (dyL+R i) + (3 y 0u;, N\) , t > 0. 


(19) 


At last, by taking the inner product of (14) with y0 u; (y) on both sides, we can get 

D ({y&u,, dyRi) - 1||0 W |||2) =w (0 W - i?i, yd^cpu) - ^ (y0 w , yi? 2 ) - £ (y/>w, #2) 

- (y<^, i'--R 2 ) + (y&>, at 2 ) , t> 0 . 

Finally, combing ( [l6| ), ( pT| ), ( p~9] ) and p0| ), defining 

cr(t) := (cj(t), v(t), D(t)^(t)) T G M 2d+2 , t > 0, 

a matrix A(t) := 


( 20 ) 


(^ 0a; , <9 cj 


< 0o;,yi?2 > 


- < d 2 (j) UJ ,R2 > | < do;0a;,y(ill + 0a;) > 


■'C 0a; 5 dyR\ ■'C 0a;? 9 ? 2 ^ 

^ dca 0a; ? dyR,2 ^ ^ ch; 0ca, 9 ?i —l - 0cu 


dujtfiuji @yR2 ^ 4 (2 < d y <f>u)jyRi ^ 11 11 ^2 ) ^ d y <f>u>i d y R 2 ^ ^ d y (/)aj>Ri ^ 

\< -Ri - 0 w ,y9a;0w) ^<y4>u,yR2> < y(j> M ,d y Ri > -\\\4> U \\ 2 L 2 <y4>u,,R2> 


for f > 0 and a column vector 


F(t) := 


-<(/) u) ,N 2 > \ 

^ tlo; 4*uj ? -^1 

< dytpu, Ni — L + R\ > 

\ < y0o;, A^2 - L-R 2 > ) 

5 


, t > 0, 









( 21 a) 


then together with ( fl4| ), we get the full modulation equations as the following coupled system 
A(t)&(t) = F(£), t > 0, 

d t R( y, t) = JH(t)R(y, t) + (£>- dy) R( y, t) + ^ JR{ y, t) - JN( y, t) 

+ J ^(y) + (-D • dy) <Uy) - w<9^(y), y € M d , t > 0 , ( 21 b) 

<t( 0 ) = (w 0 , v 0 , D 0 , £o), R( y, 0 ) = i?o(y), y e R d , ( 21 c) 

with 

N = (3 (14, + i?| 2 ) + i?) - P (<t>l) 4> u -(3 (<ti) R - 2(3' (4>l) 4>l ^ . 


In the modulation equations, ( |21a| ) is a (2d + 2)-dimensional ODE system and (21b) is a 2-dimensional 
dispersive PDE system. They are nonlinear and coupled. The components of cr which provide the information 
of the soliton, i.e. the shape, velocity, frequency and shift, are sometimes referred as collective coordinates 
in physics [39i. Note that the last component £ in a does not directly provide the original frequency 7 of 
the soliton in (A 2 ). Its derivative £ instead of itself is truly involved in the modulation equations ( | 2 l| ) and 
indicates the dynamics of 7 by ( 11 ), i.e. 


7 (t) = t(t) - - 


W 


v (s)ds - h(i) • D(t), 


t > 0 . 


( 22 ) 


Thus to give the complete information of the soliton, besides solving the modulation equations for cr, one 
also needs to solve the ODE ( 22 ) at the same time for 7. In the modulation equations ( 21 ), everything is 
real-valued. After solving it, we restore the complex-valued scalar R = Ri + iR 2 to give the radiation. 


2.2. Brief review and discussion 


The modulation equations (21) in asymptotic orders, have been studied mathematically in mm- It has 
been pointed out that when the initial radiation Ro is small, under certain conditions on the nonlinearity 
/3, the modulations equations (21) are well-posed. Moreover, the solution of the modulations equations 
preserve the orthogonality (15), and via (A), it gives the solution of the NLS equation 0 - When t 00 , 
the radiation R will vanish at the same rate as the solution to the linear Schrodinger equation with constant 
coefficients, and the collective coordinates a will turn to reach a steady state. 

We remark the modulation equations ( 21 ) are for describing the multichannel dynamics of a moving 
soliton and radiation §. For the case of a standing wave and radiation, the modulation equations are 
different and have been studied in |35] [39] from both theoretical and numerical points of view. For the 
multiple solitons with mutual interaction and radiation case, the modulation equations and the mathematical 
analysis will become very complicated, so it has not been studied a lot yet. We refer the readers to [30] for 
some results. 

The method of modulation equations described above is closely related to the method of collective 
coordinates. In fact, it was proven in [35] that if the solution at any time is close to a soliton, then it can be 
written as a soliton plus small remainder which is also orthogonal in the above sense. Hence, both methods 
give a decomposition with small corrections. However, the method of modulation equations also give a PDE 
for the radiation term, and coupling between the radiation and the ODE’s. Hence the modulation equations 
approach allows : (i) Control of the error in the ODE’s. (ii) Allows approximating the effect of radiation on 
the soliton dynamics, either exactly, or by using a good approximation of the coupling term, (iii) When the 
radiation effect is critical, as in dissipation mediated processes by the radiation [3B 1 , j37j, one can derive the 
leading dissipation (radiation mediated!) term from the couples equations, and find the leading behavior 
for processes in which a soliton changes state, for example( from excited to ground state, see [37] [38]5. 
(iv) Resolving the soliton part as it arrives the boundary of the domain of computation. This can not 
be handled by absorbing boundaries [34]. (v) Allows the rigorous asymptotic stability and scattering over 
arbitrary large time intervals. 
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3. Numerical method 


In this section, we present the numerical methods for solving the modulation equations. To do that, we 
first write down the numerical algorithm for solving the eigenvalue problem |8| for the soliton, and then we 
give the numerical method for the modulation equations (21). 

3.1. For the eigenvalue problem 

Since the modulation equations (21) replies on the solution of the eigenvalue problem (|8| for all times, 
so we need to call for a numerical algorithm for firstly finding the nonlinear bound state <fi u from (|8| 

A<A w (x) + /3 (</£(x))<A w (x) = w<^(x), x e K d . 

For a given uj > 0, we remark that the above equation is a nonlinear elliptic problem which can be solved 
numerically by the classical Newton method. However the requirement by the Newton method of an accurate 
enough initial guess will cause low efficiency in discretizing the modulation equations (21) later [39] . This 
numerical burden becomes severe especially in high dimensions. Besides, more erratic failures of the Newton 
method for solving the travelling waves have been pointed in [9j. Thus, here we apply either the Petviashvili’s 
iteration method m or the numerical algorithm which is proposed in [39] in spirit of the generalized 
Petviashvili’s method DHH, to solve ©• 

In details, for the pure power nonlinearity case in (J8|, i.e. /3(p) = A p m for p G M with A G M, m > 0, we 
apply the standard Petviashvili’s iteration method, which reads as the following. Denote ^ (n = 0,1,...,) 
as the approximation to <f> u and the Fourier transform 


02 ( k ) = [ C(x)e- ikx dx, /»(k) = [ / n (x)e _ik ' x dx, 

JR d J R d 

where / n (x) := A(</>™(x)) 2m+1 . Suppose ^ is the initial guess, then for n > 0, 

/ n (k) 


k G 


02 +1 ( k) = (M n ) c 


uj + |k| : 


k G 


(23) 


with 


M n = 


f R d G + |k| 2 ) [02( k )] 2 dk 
fud 02( k )/"(k)dk 


(24) 


and a = 2 ^ 1 which gives the fastest convergence rate of the iteration [28]. The iteration is stopped when 


\M n — 1| <e, 


(25) 


for some chosen threshold e > 0. 

The sequence {02}n>o from the above Petviashvili’s method has been proved rigorously to converge to 
the soliton (j) u in one and two space dimensions in [25]. From the application point of view, the Petviashvili’s 
method has been recognized as the most efficient numerical scheme for computing the solitary waves in a 
class of problems like with power nonlinearity and without external potential m- An improved version 
is considered in m, but it is more involved. 

For the practical implementation of the iterative method, since the bound state decays very fast to zero 
at far field [12], we truncate the problem (J8| onto a finite interval D = [— L,L] d and impose the periodic 
boundary condition for numerical issues, i.e. 

A0 w (x) + ^(^( x ))< ? !) UJ (x) = w<£ w (x), xed, (26a) 

<Mxj) = 0w(x_j), e dfl, (26b) 
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where x = (x lt .. .,x d ) T , xj = (x lt ... ,Xj-i,L, x j+1 ,. ..,x d ) and x_j = (xi,.. .,xj- i, -L,x j+1> . ..,x d ) for 
j = 1,..., d. The Fourier transforms and integrations in the Petviashvili’s method (23)-(24) are implemented 
by means of the Fourier pseudo-spectral method [42] . For example, in one dimension Q = [—L,L], if we 
denote Xj = — L + j • ^ for j = 0,1,..., N with N an even integer, then 

or N — l N — 1 7 

#;(*o*^ E rjxj)e- ik ^ +L \ nh) (2?) 


i=o 


j=o 


for l = —iV/2, -JV/2 + 1,..., 7V/2 - 1, and 


M 77 


Eri -^/2 ("+*?) i ^(^)] 2 

-\JV/2—1 


(28) 


EJ-w/ 2 0S(*O/ n (*O 

This implementation introduces some extra error. The choice of the L determines the boundary truncation 


error of the soliton and the error of the quadrature for approximating the integrals as in (28). The choice 


of the N controls the interpolation error of the Fourier transforms as in (27). Both error should decay 
exponentially when L and N increase, if the solution is smooth and decay very fast at the far field. 

For general nonlinearity case in ([8]), i.e. /?(■): M —)> M is some general nonlinear function, we apply the 
numerical algorithm proposed in [39|. By adopting the same notations as above and truncating the problem 


as (26), the algorithm reads: 


Step 1 Find the ground state of the Hamiltonian functional 

H n { 4 >)= [ [|V<+-/ 3 ((C) 2 ) M 2 ]dx, 

Jn 

in the unit sphere of L 2 (D). Denote the solution as 

C +1 := axgmm{H n ((J)) : 0 e L 2 (fl), \\cf >\\ L 2 = 1, +x) > 0}. 


(29) 


Step 2 Scale the ground state </>2 +1 according to the energy uj. That is to find the scaling constant c n 
such that 

C +1 := c n C +1 , (30) 

satisfying 


[ [-m 

Jn 


n+l 1 2 


+ ((C +1 ) 2 ) (C +1 ) 2 ] dx = w ||c +1 ||^, 


which is obtained by taking the inner product of (26a) on both sides with in L 2 (S4). Then we can solve 
the equation 


Jj ((c"c +1 ) 2 ) ^- +1 ) 2dx=UJ+ i 


v# 


n+l 


dx, 


for the value of c n . In particular, when it comes back to the power nonlinearity case /3(p) = Ap m , A G M, m > 
0, we have explicit formula 


W + Jfi 

vr E +1 

2 

dx 

^ In 

r E +1 

2m+2 

dx 


Then iterate until o converges. The Cauchy criterion is used as the stopping condition, i.e. 

IIC +1 -CllL~<e, (31) 

with some chosen threshold e > 0. 
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To implement the above algorithm, for the first step, we can use the normalized gradient^ flow method 
with a backward Euler Fourier pseudospectral discretization m [7] to get the ground state . For the 
second step, we use the standard Fourier pseudospectral discretization [42] for the spatial derivative and 
integrations. Once is obtained, d can be found out from (10) with periodic boundary conditions on 
It by the Fourier pseudospectral discretization again, and so doesdfcfu in (18). The rigorous mathematical 
analysis of the convergence of the above algorithm (29)-(30) has been given in [39]. The algorithm has been 
shown to be efficient, robust even if the initial guessl/r, is quite away from the exact solution, and converge 
to the (j) u rapidly. We refer the readers to [39] for more details. 


Remark 3.1. There are many generalized versions of the Petviashvili’s method proposed in the literature 
in order to deal with the eigenvalue problems involving different kinds of external potentials, the coupled 
equations case and the general nonlinearity case J7J fT/i I ~23f . Most of them are proposed with the ad-hoc 
approximation. The scheme varies when it comes to different cases and it would need choices of some 
free parameters. They are shown numerically to converge well for some specific problems, but the general 
convergence results of them are not available. While, the algorithm proposed in does not have 

the ad-hoc regularization. The normalization is done in a precise way each step. The scheme is defined very 
precisely for general cases and the convergence is guaranteed mathematically in 139]. This is the reason why 
we use the algorithm rather than the others such as the spectral renormalization method m ■ 


Remark 3.2. The iteration algorithm (29)-(30) also applies to the pure power nonlinearity case in \2b 
It appears to more complicated than the Petviashvili’s method (2tf\)-(24), since it is originally derived in 
m with an external potential where the linear differential operator and the potential have to be treated 
separately in the computation. In that case, it is impossible to stay just in the Fourier frequency space for 
computing. Another issue that makes (29)-(3Cfy more involved than the p^-p^]) is the call of an imaginary- 


time evolution method for obtaining the ground state of linearized Hamiltonian. Although we are using the 
classical normalized gradient flow method for simplicity, one needs to realize there are many new developed 
numerical techniques in recent research for rapidly increasing the efficiency of computing the ground state. 
Thus, a fair comparison between the two iterative methods would need a systematical study, which is beyond 
the scope of this paper. 


3.2. For the modulation equations 

Similar to the eigenvalue problem, si nce t he dispersive wave R also decays very fast to zero at far field 
at finite time, we truncate the problem (21b) onto a finite interval ft = [—L,L] d and impose the periodic 
boundary condition for computation. The truncated initial boundary value problem of the modulation 
equations read, 


A(t)a(t) = F(t), t > 0, (32a) 

d t R( y, t) = JH(t)R(y, t) + (i) • dy) R(y, t)+(j+ JR( y, t) - JN( y, t) 

+ (f + J K(y) + [d • d y ) K{y) - udu&ui y), y e n, t> o, (32b) 

cr(0) = (cj 0 , v 0 , D 0 , Co), R( y,0) = -Ro(y), y e ft, (32c) 

R(yj)=R(y-j), d Vj R{y,) = d yj R{ y_ j), yj,y-j edD, j = l,...,d, (32d) 

where yj,y_j are defined similar as Xj,x_j. 

Choose the time step size r = At > 0 and denote the time steps by t n := nr, n = 0,1,.... To present 
the scheme, we denote 


a n = (cu n ,v n ,D n ,C)~a(t n ), A n « A(t n ), F n « F(t n ), 
H n ^H{t n ), 7 n «7 (tn), R n {y) ~ R(y,t n ), 

C(y) ~ ^w(t„)(y)> d u 4>2(y) & d u <l> u ( tn )(y), N n (y) « N(y,t n ), 
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and introduce the finite difference operator on some grid functions / n , 

jn -\-1 _ yn—1 


s t f n ■■= 


2 T 


n = 1 , 2 ,.... 


Then a semi-implicit leap-frog finite difference temporal discretization of (32) reads, 

A n S t a n = F n , n = 1 , 2 ,..., 

StR n ( y) = ^ (Jff n + <5 t £> n • d y ) (f? n+1 (y) + R n ~\ y)) + (0C + • y^ JfC(y) 


(33a) 


with initial values 


- JN n (y) + + -<5 t v n • yj JC(y) + (6 t D n • d y ) C(y) 

- <5 t w n <9 w </>”(y), yen, n=l,2,..., 


a 0 = (cj 0 , v 0 , -D 0 , Co), R°(y) = Ro(y)- 


(33b) 


Also, together with (33) we dicretize (|22| as 

n 

^ 7 " = StC -jYl S ^ n ■ - o 6 ^ n - Dn ’ n = °, 1 , - -, 


m —0 


where we apply the composite trapezoidal rule to approximate the integral term. Since (33) is a two-level 


scheme, we also need the starting values at t = t\. In order to get a second order accuracy in temporal 
approximations, they are obtained by the Taylor’s expansion of the solution and noticing the equations 
(33b|) as 


a 1 = a 0 + rS t o-° , S t cr 0 = (A 0 ) 1 F°, 

7 1 = 7° + rSf 7°, -R^y) = -R°(y) + T5 t R°(y), 

S t R°( y) = ( JH° + 6 t D° • d y ) R°( y) + + ^ t v° • y^ JR°(y) 

- JN°(y) + Ac 0 + ^<5 t v° • y) (y) + (<5 t D° • 9 y ) 0° (y) - (y). 


(33) is the semi-discretization of (32). To get the full discretization, i.e. to discretize the space and ap¬ 
proximate the above spatial derivatives, we use the standard Fourier pseudospectral method [42]. Thus, our 
numerical method can be referred as the semi-implicit Fourier pseudospectral (SIFP) method. Here, 02 is 
obtained by algorithm (29)-(30) from (26) and 9^02 i s gi ven by (10). 

The SIFP method is clearly time symmetric. In the scheme of SIFP, (|33a|) is fully explicit, while (33b) 


is semi-implicit. So at each time level t = £ n , we apply a linear solver, for example the Gauss-Seidel method 
PI, to get R n+1 . We re mark that here the reason why we put a time average on the function behind the 
operator JH + D • d y in (33b) is to get rid of the stability problems [6]. Finally, we would like to comment 


that this modulation equations approach as a numerical solver for the NLS equation works in any space 
dimensions, since the modulation equations are consistent with the NLS and the corresponding numerical 
discratization are well-defined in any dimensions. 

Remark 3.3. We remark that in the algorithm for the eigenvalue problem (26) in Section [3] and 

the temporally discretized modulation equations \3cfy, one can also use the finite difference method for spatial 
discretizations. Here we choose the Fourier pseudospectral method for a high accuracy purpose in the case 


of periodic boundary conditions (26b) and (3 2d). 

Remark 3.4. In the imposed boundary condition is the zero boundary where corresponding sine spectral 
method is applied, while here we use the periodic boundary condition (26b) and \32df). Both types of boundary 


conditions are fine to use as approximations to the physical model after domain truncation if the finite domain 
is chosen large enough. The reason why we use the periodic boundary condition here is because we still want 
to apply the spectral method in the presence of the gradient operator in ( 32b\). 
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4. Numerical results 


In this section, we present the numerical results of modulation equations via using the proposed numerical 
method. To do that, we first test the correctness and accuracy of the SIFP method. Then we apply the 
SIFP method to numerically explore the multichannel dynamics. 


4-1. Accuracy test 

For simplicity, we consider the one-dimensional case, i.e. d 


1 and 


x = x, 


V = V, y = y, 


in the NLS equation 0 and the modulation equations ( |2l| ), to test the SIFP method ( [33] ) . We take the 
cubic nonlinearity, i.e. 

/3(\u\ 2 )u = X\u\ 2 u, A = 2, (34) 

in 0 and choose the computation domain ft = [—L, L\ for the variable y in ( [32] ) with L = 16 which is 
large enough to ignore the boundary truncation error before the wave reaches the boundary of the ^/-domain 
during a short time computing. We choose the initial data in 0 

vo = 0.5, ujo = 1 , 70 = 1 , D 0 = 0 , 

where we know explicitly the corresponding (j) UQ = sech(x). Then in order to satisfy the orthogonality 


conditions (15), we choose R{pc, 0) = Ro(x) = Ri(x,0) + iR 2 (x,Q) as 

< 0.2xe~ x2 > 


Ri(x, 0) = 0.2xe 


II x<p, 




R 2 (x,0) = 0. 


(35) 


WQ II L 2 


The threshold 6 used for the stopping criterion ( [ 25 ] ) for the iteration algorithm (29])- (30) is chosen as 5 = 10 -9 . 

To show the modulation equations with SIFP solve the NLS equation jl) correctly, we solve the mod¬ 
ulation equations (32) numerically by the SIFP (33) to get 0 w m(i/), v M , y M , D M and c o n for 

0 < n < M = T/r, and use the ansatz (A) to construct the numerical solution u M (x) of the NLS equation 
0 with application of the composite trapezoidal rule to approximate the integrals, i.e. 


vl (x) :=exp i 



Then, we compute the error 

e u (x, T) := u(x, T) — 4> M (x), xeft, (36) 

where the exact solution u(x, T ) of the NLS equation 0 is obtained by classical numerical methods such as 
the time-splitting Fourier spectral method mil] with very small step size, e.g. r = 10 4 , h = 1/16. We test 
the temporal and spatial discretization errors of the SIFP method separately. Firstly, for the discretization 
error in time, we take a fine mesh size h — 1/8 such that the error from the discretization in space is negligible 
compared to the temporal discretization error. The errors (36) under maximum norm are presented at T = 1 
in Tab. [l] Secondly, for the discretization error in space, we take a very small time step r = 10 -4 such 
that the error from the discretization in time is negligible compared to the spatial discretization error. The 
corresponding errors under maximum norm are presented at T = 1 as well and tabulated in Tab. [2] 

Clearly, from Tabs. [l][2] we can conclude that the SIFP method (33) solves the multichannel solutions 
based on the modulation equations for the NLS equation 0 correctly and accurately. The numerical method 
has second order accuracy in time and spectral accuracy in space. 
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Table 1: The temporal error and convergence rate of the SIFP method for the modulation equation under 
different time step r with h = l/8atT = l. 



t —1 

0 

II 

P 

7 - 0/2 

to/2 2 

to/2 3 

ro/2 A 

ro/2 5 

\\ e u II Loo 
rate 

8.70E-3 

3.20E-3 

1.44 

1.20E-3 

1.46 

3.78E-4 

1.66 

9.76E-5 

1.96 

2.41E-5 

2.01 


Table 2: The spatial error of the SIFP method for the modulation equation under different mesh size h with 
r = 10- 5 at T = 1. 



h 0 = 1 

h 0 /2 

ho/2 2 

h 0 /2 3 

\\ e u\\ l°° 

2.40E-2 

8.03E-5 

4.02E-8 

2.55E-9 


4.2. Comparisons 

Now we compare the modulation equations approach to solve the dynamics in the NLS equation with 
existing direct numerical methods toward discretizing q. We use the same setup as in ( |4.1| ) but with a 
larger velocity vq = 16. We work on the fixed computational domain Q = [—16,16] by using the SIFP 
method for the modulation equations (32) with the translated variable y. The dynamics of the multichannel 
solution in the original x-domain are shown in Fig. |T] In Fig. [ 2 J the same NLS equation ([!]) problem is 
solved directly by imposing zero boundary condition and using the time-splitting sine spectral method EJ[8]. 

Based on Figs. [l]and[2j we can see that on a domain of the same size, within the computational time, 
the radiation wave in the modulation equations has not reached the boundary yet in the ^-domain, but the 
solution of the NLS equation has already hit the boundary due to the velocity. The waves in the modulation 
equations approach can pass through the boundary of the x-domain naturally, while the waves in the NLS 
equation are destroyed by the zero boundary condition and the direct PDE solver on x-domain. We remark 
that even if the radiation wave reaches the boundary of ^/-domain, we can use the absorbing boundary 
techniques to the equation (32b) of R and improve the results. 

Of course for this simple case, one can change to the periodic boundary condition to avoid the breakdown 
of the wave at the boundary. The following example illustrate the fails of the periodic boundary condition. 

When there are two solitons in the multichannel dynamics and two solitons are well-separate initially 
and moving towards opposite directions, the interactions between the solitons can be ignored. That is to 
say they can be treated as two single solitons in the NLS equation which can be handled by the modulation 
equations (21). Here we choose the same numerical example but with the velocity and shift 


^0 = ±7, Dq = ±7, 


then the dynamics of the multichannel solution is shown in Fig. [3j We also show the corresponding results 
of solving the NLS equation directly by imposing periodic boundary condition and using the time-splitting 
Fourier spectral method. 

By this numerical experiment, based on the results in Fig. [3j we can see that the modulation equations 
approach can simulate the multiple solitons case if the solitons are always well separated during the dynamics. 
However the direct solver for the NLS equation with periodic boundary condition can not, because the waves 
will collapse when they hit the boundary and enter the domain from the other side. 

Now we apply the SIFP method to study the dynamics of the multichannel solutions with the same setup 
as in (4.1) numerically. In order to provide a long time simulation, we enlarge the domain to = [—64, 64], 
such that the dispersive wave R is always away from the boundary during the simulation. Take step size 
r = IE — 3, h = 1/8, we solve the modulation equation ( [32] ) till the collective coordinates cr(t) reach the 
steady state. The dynamics of the collective coordinates v(t ), D(t ), uo(t) and y(t) are shown in Fig. [4] The 
profiles of the dispersive wave R(y,t) at different time are shown in Fig. [ 5 ] 
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Figure 1: Dynamics of the solution \u(x,t)\ to the NLS equation in x-domain and the <j) u (y) and R(y,t ) at 
t = 1 in the ^/-domain. 



o L_-*-— —-*— - -— 

— IS 1 O —5 O 5 lO 15 
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Figure 2: The solution \u(x,t)\ to the NLS equation in the x-domain Q = [—16,16] by using zero boundary 
condition and the time-splitting spectral method. 


4-3. Explorations on multichannel dynamics 

With Dq = —25 the long time dynamics of the multichannel solution in the numerical example (34)-(35) 
on a large domain Q = [—1150,1150] is shown in Fig. [6j 

Based on Figs. |4pl we can draw the following observations: 


1. The collective coordinates v(t), D(t ), cj(t) and y(t) converge the steady state as t goes to infinity (cf. 
Fig. [4|, and the dispersive part R(y,t) spreads out to far field (cf. Fig. [5|. 
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Figure 3: Dynamics of the solution | u(x, t)\ to the NLS equation in x-domain: The results of the modulation 
equations method (left figure); The results of time-splitting spectral method (right figure). 





Figure 4: Dynamics of the collective coordinates cr(t). 


2 . 

3. 


The dynamics of each component of <j(t) is not monotone in time t (cf. the left figure in Fig. [§. This 
indicates that the process of the dynamics of the soliton and the dispersive wave is not monotone. 
The dispersive wave 0 in this case has a large expanding velocity (cf. Fig[5|. It is because that the 
chosen initial perturbation, i.e. the Rq in (35), has a large H 1 norm. For the kind of situation, we 
remark that using the absorbing boundary conditions could be a more efficient way of study. 


14 








































Figure 5: Profiles of the dispersive wave \R(y,t)\ at different t. 


The modulation equations method also works well in two dimensions. Here we give a 2D numerical 
example. We take d = 2 and 


x = (xi,x 2 ), y = {yi, y 2 ), v = (ui,i; 2 ), D = (D 1 } D 2 ), 


m 


0 and ( [2l] ) with cubic nonlinearity 


We choose the initial data as 


Vo = ( 1 , 0 ), D 0 = (0,0), 70 = 0.5, w 0 = 1, 


and 


i?i(x,0) = 0 . 2 xie“^"^ 


< 0.2xie X *,xi 4>u 0 > 

\\Xl<t>LO 0 \\ 2 L 2 




i? 2 (x, 0) = 0. 


and the computational domain as H = [—16,16] x [—16,16]. The dynamics of the soliton 0 w (x, t) : 
4> u (t)( y(x, t),t) and the radiation i2(x, t) := |i2(y(x, t), t)| by the modulation equations are shown in Fig. 
with camera fixed in domain [—8,8] x [—8,8]. 
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